Electronic stopping power in gold: The role of d electrons and the H/He anomaly 
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The electronic stopping power of H and He moving through gold is obtained to high accuracy 
using time-evolving density-functional theory, thereby bringing usual first-principles accuracies into 
this kind of strongly coupled, continuum non-adiabatic processes in condensed matter. The two key 
unexplained features of what observed experimentally have been reproduced and understood: (i) 
The non-linear behaviour of stopping power versus velocity is a gradual crossover as excitations tail 
into the d-electron spectrum; and (ii) the low- velocity H/He anomaly (the relative stopping powers 
are contrary to established theory) is explained by the substantial involvement of the d electrons in 
the screening of the projectile even at the lowest velocities where the energy loss is generated by 
s-like electron-hole pair formation only. 



Non-adiabatic processes are at the heart of aspects of 
science and technology as important as radiation dam- 
age of materials in the nuclear and space industries, 
and radiotherapy in medicine. Yet, in spite of a long 
history, the quantitative understanding of non-adiabatic 
processes in condensed matter and our ability to per- 
form predictive theoretical simulations of processes cou- 
pling many adiabatic energy surfaces is very much behind 
what accomplished for adiabatic situations, for which 
first-principles calculations provide predictions of var- 
ied properties within a few percent accuracy. Substan- 
tial progress has been made for weakly non-adiabatic 
problems such as the chemistry of vibrationally excited 
molecules landing on metal surfaces [1], but not in the 
stronger coupling regime of radiation damage. Recently, 
the electronic stopping power for swift ions in gold has 
been carefully characterized by experiments [2H1], show- 
ing flagrant discrepancies with the established paradigm 
for such problems [5, 6 , and only qualitative agreement 
with time-dependent tight-binding studies [7 , and with 
detailed studies for protons based on first principles [5], 
leaving very fundamental questions unanswered in spite 
of the apparent simplicity of the system. Most notably 
the H/He anomaly: the present understanding predicts a 
stopping power for H higher than for He at low velocities 
[5] , which strongly contradicts the recent experiments [I] . 

A particle moving through a solid material interacts 
with it and loses its kinetic energy to both the nuclei 
and the electrons inside it. At projectile velocities be- 
tween 0.1 and 1 atomic units (a.u. henceforth) both the 
nuclear and the electronic contributions to the stopping 
power (energy lost by the projectile per unit length) are 
sizeable [7|. Based on the jellium model (homogeneous 
electron gas) the electronic stopping power, S e , is pre- 
dicted to be S e oc v for a slow projectile traversing a 



metallic medium [HI [TU] . Such behaviour has been ob- 
served experimentally in many sp-bonded metals |lllll2j . 
and the jellium model has allowed deep understanding of 
the dynamic screening of the projectile and its relation to 
stopping |13j . Even the jellium prediction of an oscilla- 
tion of the proportionality coefficient with the projectile's 
atomic number Z has been verified Jfij and reproduced by 
ab initio atomistic simulations |14| . However, phenomena 
that cannot be accounted for within the jellium paradigm 
have been described only qualitatively so far [TSHTS] . Ex- 
periments on noble metals Cu, Ag and Au, show pro- 
nounced nonlinearities in S e (v) [51 ITT1 IT51 HTJ] . In the 
case of slow H and He ions in gold [2H1], S e (v) displays an 
increase in the slope roughly around v ^ 0.18 a.u. This 
is usually attributed to a threshold projectile velocity 
needed to excite the d-band electrons that are relatively 
tightly bound. A model was developed based on the ab 
initio density of electronic states and a stochastic treat- 
ment of excitations |17j . which reproduces the threshold 
for protons. Here we obtain the non-linear S e (v) and the 
H/He anomaly with a general purpose ab initio method 
equally applicable to many other radiation problems. 

We calculate the uptake of energy by the electrons in 
gold from a moving H or He ion in its (100) channel by 
explicitly following the dynamics of the electrons cou- 
pled to the projectile's motion, using time-evolving time- 
dependent density functional theory (TDDFT) [20]. We 
find very good quantitative agreement with some recent 
low energy ion scattering experiments on thin gold films 
[5H3] . The results are analyzed in terms of the electronic 
excitations that are responsible for the energy loss, which 
very clearly shows why the slope of S e increases with pro- 
jectile velocity. In contrast to the usual idea that at low 
projectile velocity only electrons close to the Fermi en- 
ergy contribute to the stopping, we find that there is a 
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significant contribution from deep lying states even for a 
slow projectile. This means that at low velocities (v < 0.2 
a.u.) the electrons accessible to excitations (s) are differ- 
ent from the ones involved in the screening of the projec- 
tile (s+d), the latter providing the excitation mechanism. 

We performed all calculations using the Siesta 
method [21] in its time evolving TDDFT implementa- 
tion [55]. We used the Perdew-Burke-Ernzerhof (PBE) 
version [33] of the Generalised Gradient Approximation 
(GGA) to the instantaneous exchange and correlation 
functional. Since we are only interested in the very low 
energy regime that is far below core electrons excitation 
thresholds, only valence electrons in Au are considered 
explicitly and a norm-conserving pseudopotential is used 
to describe the core electrons (up to the 5p sub-shell). 
Further details are found in [24 . After placing the pro- 
jectile in a [001] channel and finding the DFT ground 
state, it is given an initial velocity along the z-direction 
while all gold atoms are initially quiescent. The system 
evolves by following Ehrenfest coupled electron-ion dy- 
namics. On the time scale of the simulation (~ 0.75 — 6.0 
fs for v = 0.05 — 0.50 a.u.), the gold nuclei only gained 
negligible velocities and did not move significantly. We 
monitored the total energy of the electronic subsystem 
as a function of time. Once the transient related to the 
sudden start has disappeared [T51 [M] , S e is extracted as 
the average rate of change of the electronic energy with 
the distance travelled by the projectile. 

Fig. [I] shows our results for S e (v) for H and He projec- 
tiles in gold for the velocity range v — 0.06 — 0.50 a.u. We 
also plot results of some recent experiments performed 
on thin single crystal gold films oriented along (100) [2] 
and polycrystalline gold films [3] [4] . The agreement be- 
tween our simulations and the experiments is noticeable. 
Although the stopping power is still underestimated (es- 
pecially for H around v — 0.3 a.u.), no previous ab initio 
approach had this level of agreement on the non-linear 
velocity dependence of the stopping power of real mate- 
rials. Our results for the stopping power are well con- 
verged with respect to the basis size and the density of 
points on the real and the momentum space grids |24j . 
A larger basis set for the projectile, however, (TZDP in- 
stead of DZP 24J) increases the stopping power about 
5% at v = 0.5 a.u., but considerably less at low velocity. 
The error bars in Fig. [T] indicate the dispersion in our 
results for v =0.08, 0.1 and 0.5 a.u. when the various pa- 
rameters are varied, including the basis set [21] (the bars 
for low velocities are hardly larger than the size of the 
circles). The strict channelling in the simulation is partly 
behind the observed underestimation: calculations for a 
30% smaller impact parameter give a 25% increase in S^ 1 
at v = 0.28 a.u. that reduces to 1% at v = 0.5 a.u. 

We see a clear deviation from the linear behaviour 
around v — 0.2 a.u. in S e of both H and He. This is 
unlike the S e cx v of the uniform electron gas. It seems a 
plausible explanation that at low projectile velocity only 
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FIG. 1. Electronic stopping power of H and He projectiles 
in gold as a function of projectile velocity. Results of our 
simulations are compared with the experimental data from 
Refs. [2H1] on single and polycrystalline thin gold films. 



s-band electrons from the states around the Fermi energy 
contribute to the stopping and at higher velocity elec- 
trons in the d band that lie relatively deeper in energy 
are also able to take part in it, resulting in an increase 
in the slope of S e . Thus, comparisons have been made 
[2 0] with jellium using the average electron density n e 
of the s electrons (r s — 3.01 a.u., where n^ 1 — §7rr;?), 
using r s — 1.49 a.u., corresponding to an effective num- 
ber of s and d electrons [3], or of the density in the (100) 
channel (r s = 1.8 a.u.). However, the jellium predictions 
do not agree with the experimental results except at pro- 
jectile velocities around v = 0.6 a.u. in the latter case, 
despite the expectation that all the d-b&nd electrons are 
active for a projectile velocity v > 0.47 a.u. [TS]. There 
is a further problem in the comparison with jellium: If 
we assume that at low velocity only s electrons are ac- 
tively participating in the stopping mechanism, the jel- 
lium model predicts Sf > Sf e [5], which is not the case. 

To explain the above inconsistencies and get a bet- 
ter idea of the energy loss mechanism we compute the 
changes in the electronic distribution due to the ex- 
citation of the electrons when a projectile propagates 
through the material. Having {\ip n {t))} and X(t), the 
set of evolved occupied KS states, and the correspond- 
ing atomic positions at time t, we calculate the adiabatic 
states {\<f>i,X)}, i.e., the set of self-consistent static KS 
states for X(<). By projecting the evolved states onto 
the adiabatic states, C„ = (cf)i,'K(t)\ip n (t)) , the density 
of occupied energy states 0(E) at time t as a function of 
energy E are obtained as 0(E) — J2in \Cin\ 2 8(E — ^»)- 
Here Ei is the eigenvalue of the adiabatic state |0j,X). 
To compute the change in the electronic distribution or 
the (electron-hole) excitation distribution, P(E), we sub- 
tract the ground state electronic distribution from 0(E). 
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FIG. 2. (Color online) Excitation distribution P(E) when 
H passes through gold with velocity 0.05 a.u., for time values 
between t = 0.1 fs and 1.1 fs in steps of At = 0.1 fs (light color; 
larger amplitude for longer t; Gaussian broadening a — 0.2 
eV). The dark curve is the electronic density of states g(E) 
(a — 0.5 eV). P{E) and g{E) are in different scales. 



That is, P(E) = 0(E)-Q(E F -E)g(E), where E F is the 
Fermi energy of the system, g(E) is the electronic density 
of states and Q(E) is the Heaviside step function. 

Fig.[2]shows the excitation distribution P{E) as a func- 
tion of energy at various instants from t = 0.1 fs to t = 1.1 
fs with an interval At = 0.1 fs, for the passage of a H 
atom in gold along (100) with velocity v = 0.05 a.u. The 
electronic density of states of the bulk Au host g(E) is 
also plotted in Fig. [2j The negative and the positive val- 
ues of P{E) show the density of empty and filled states 
below and above Ep , respectively, due to the electronic 
excitations caused by the moving projectile. Notice that 
despite being very slow (v — 0.05 a.u), the projectile is 
able to excite relatively tightly bound <i-band electrons. 
A short initial transient behaviour is noticeable in Fig. [2] 
at energies deep below the Fermi energy, the number of 
empty states becomes larger initially, requiring a short 
time to adjust to a stationary regime. This is because in 
our simulations the projectile is a static impurity atom 
at t = that suddenly acquires a finite velocity resulting 
in a large initial perturbation. 

To see how the excitation distribution after the tran- 
sient depends on the velocity of the projectile, we plot 
P(E) against E in Fig. [3] at t = 0.25 fs for various pro- 
jectile velocities, between 0.05 a.u. and 0.50 a.u. We see 
that, compared to the states just below the Fermi energy, 
the number of excitations from deep inside the d-band in- 
creases more quickly with the velocity of the projectile. 
This means that the effective number of d-band electrons 
involved directly in excitations provoking the stopping 
process increases with the projectile velocity. To see this 
more clearly, we separated the energy window into two 
parts at the upper edge of the d-band at energy Ed and 



calculated the total number of excitations N% and N 2 
from the states below and above Ed for a constant dis- 
tance travelled by the projectile. We find that P(E) oc t 
after the initial transient so we can estimate Ni and N 2 as 
^1 = 5 I-L \P{E)\dE and N 2 = A f£ \P(E)\dE, which 
we did using P(E) at t = 0.25 fs. In the inset of Fig.[3]we 
plot iVi and N 2 and the fraction Ni/(Ni + N 2 ) — Ni/N 
against the projectile velocity as dashed, dotted and solid 
lines for H and He projectiles. We see that N% and Ni/N 
increase with v for both projectiles. For H, increases 
and saturates whereas for He it increases up to v = 0.3 
a.u. but decreases for a faster projectile. Since there is 
one s electron and ten d electrons and Ni also includes 
the contribution from the s-band states, ideally it should 
tend to Nt/N - 10/11 = 0.909 for high projectile ve- 
locity. Ni/N reaches only 0.88 and 0.78 for H and He 
at v = 0.5 a.u. Although the fraction of excitations from 
the deep lying states is higher for H, the absolute number 
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FIG. 3. Up: Excitation distribution P(E) due to the passage 
of a H (top) or He projectile (bottom) in gold evaluated at 
t = 0.25 fs for various projectile velocities, v = 0.05 — 0.50 a.u. 
in steps of 0.05 a.u. Increased projectile velocity gives curve 
with larger amplitude (indicated by arrows). The dashed and 
dotted vertical lines show the upper edge of the gold's 5d-band 
Ed and the Fermi energy Ep. Down: Number of empty states 
below and above Ed, Ni and N2, and fraction Ni/(Ni + N2) 
versus projectile velocity, due to the excitations for H or He. 
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FIG. 4. Up: The curve shows R = Sf c /Sf for jellium 
versus the electron density parameter r 3 The values of R 
obtained for Au for v — 0.1 and 0.5 a.u. are associated with 
r a = 3.04 and 1.49 a.u, respectively, following Ref. [1]. The 
calculated ratio R for a system made of Na atoms in bulk Au 
positions, r s = 3.04, is also presented. Down: Projection of 
KS states for the system with projectile onto the KS states of 
bulk Au, for H and He, subtracting the Au density of states. 



is lower, as can be seen in the figure. Furthermore, in the 
case of H, Ni > N2 for the whole velocity range shown 
whereas for He, N2 > Ni in the very low velocity range. 

We address now the low- velocity H/He anomaly. Fig. [4] 
presents the ratio R = Sf e /Sf in jellium [5]- The values 
of R for Au for v = 0.1 and v = 0.5 a.u. are plot- 
ted on the two dotted vertical lines at r s = 3.04 and 
r s = 1.49 a.u., which correspond to 1 and 8.24 electrons 
per bulk unit cell, i.e., the s electrons and the effective 
number of valence electrons (s and d) that fit the plas- 
mon pole for bulk Au 0], We see that for the faster 
projectile R is close to the jellium value and significantly 
larger than 1. However, for the slower one, we obtain 
R = 4.7, in clear disagreement with the jellium value of 
0.79, but in agreement with experiment. We also plot R 
for the fictitious system built by putting Na atoms in the 
Au positions, which corresponds to an electron gas with 
r s = 3.04. The plot shows a perfect agreement for the R 
values of Na and jellium [25]. These differences between 
jellium (or Na) and gold are thus due to the presence 
of gold's d electrons. This is consistent with the fact 
that, even if a slow projectile were unable to excite the 
d-band electrons appreciably, the presence of the projec- 



tile in gold constitutes a large static perturbation for the 
d electrons. This can be clearly seen by calculating the 
projection of the ground state of the gold with the pro- 
jectile onto that without it and obtaining a distribution 
analogous to P(E), now describing the static screening of 
the projectile, i.e., projecting the wave-functions of Au 
with the projectile onto the states of pure Au (Fig. [2]). 
This means that for a slow projectile the response of the 
electrons in gold is far from the one described by the 
homogeneous electron gas model that includes just the 
s-band electrons. 

To summarize, we have shown that realistic non- 
adiabatic stopping of projectiles in real metals can now be 
described from first-principles with acceptable accuracy, 
even at the Ehrenfest dynamics level. We used it to cal- 
culate the electronic energy loss on passage of H and He 
through Au and find good quantitative agreement with 
experiments. Many problems can now be addressed with 
this technique in the fields mentioned in the introduction. 

JK thanks the Wellcome Trust for a Flexible Travel 
award. MAZ acknowledges support of the Islamic Devel- 
opment Bank. DSP, JIJ and AA acknowledge support 
of the University of the Basque Country UPV/EHU, un- 
der Grant IT-366-07, and of the Spanish Ministerio de 
Ciencia e Innovation under Grant FIS2010-19609-C02- 
00. The calculations were done using the CamGrid high- 
throughput facility of the University of Cambridge. 



[9 
[10 

[11 



[12 
[13 



N. Shenvi, S. Roy and J. C. Tully, and refs. therein, Sci- 
ence 326, 829 (2009). 

E. A. Figueroa, E. D. Cantero, J. C. Eckardt, G. H. 
Lantschner, J. E. Valdes, and N. R. Arista, Phys. Rev. 
A 75, 010901 (2007). 

S. N. Markin, D. Primetzhofer, S. Prusa, M. Brunmayr, 
G. Kowarik, F. Aumayr, and P. Bauer, Phys. Rev. B 
78,195122 (2008). 

S. N. Markin,D. Primetzhofer, M. Spitz, and P. Bauer, 
Phys. Rev. B 80, 205105 (2009). 

P. M. Echenique, R. M. Nieminen and R. H. Ritchie, 
Solid State Comm. 37,779 (1981). 

P. M. Echenique, F. Flores and R.H. Ritchie, Solid State 
Physics 43, eds. H. Ehrenreich and D. Turnbull, Aca- 
demic Press, p. 230 (1990). 

C. P. Race, D. R. Mason, M. W. Finnis, W. M. C. 

Foulkes, A. P. Horsfield, and A. P. Sutton, Rep. Prog. 

Phys. 73, 116501 (2010), and refs. therein. 

E. D. Cantero, G. H. Lantschner, J. C. Eckardt, and N. 

R. Arista, Phys. Rev. A 80, 032904 (2009). 

R. H. Ritchie, Phys. Rev. 114, 644 (1959). 

M. Kitagawa and Y. H. Ohtsuki, Phys. Rev. B 9, 4719 

(1974). 

J. E. Valdes and G. Martmez-Tamayo, G .H. Lantschner, 
J. C. Eckardt and N. R. Arista, Nucl. Instrum. Meth. B 
73, 313 (1993). 

G. Martmez-Tamayo, J. C. Eckardt, G. H. Lantschner, 

and N. R. Arista, Phys. Rev. A 54, 3131 (1996). 

M. Quijada, A. G. Borisov, I. Nagy, R. Dfez Muino, and 



5 



P. M. Echenique, Phys. Rev. A 75, 042902 (2007). 
[14] R. Hatcher, M. Beck, A. Tackett, and S. T. Pantelides, 

Phys. Rev. Lett. 100, 103201 (2008). 
[15] J. E. Valdes, J. C. Eckardt, G. H. Lantschner, and N. R. 

Arista, Phys. Rev. A 49, 1083 (1994). 
[16] P. Vargas, J. E. Valdes and N. R. Arista, Phys. Rev. A 

53, 1638 (1996). 
[17] J. E. Valdes , P. Vargas and N. R. Arista, Phys. Rev. A 

56, 4781 (1997). 
[18] J. M. Pruneda et al. Phys. Rev. Lett. 99, 235501 (2007). 
[19] R. Blume, W. Eckstein, and H. Verbeek, Nucl. Instrum. 

Methods 168, 57 (1980); 194, 67 (1982). 
[20] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 

(1984). 

[21] P. Ordejon, E. Artacho, and J. M. Soler, Phys. Rev. B 



53,10441 (1996); J. M. Soler et al., J. Phys. Condens. 
Matter 14, 2745 (2002). 

[22] A. Tsolakidis, D. Sanchez-Portal, and R. M. Martin, 
Phys. Rev. B 66, 235416 (2002). 

[23] J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. 
Lett. 77, 3865 (1996). 

[24] See details in supplementary materials. 

[25] The stopping power values for H and He in this Na are 
3.902 eV/A and 3.083 eV/A, respectively, for v = 0.5 
a.u., and 1.024 eV/A and 0.715 eV/A, for v = 0.1 a.u.; 
the corresponding values for jellium as extracted from 
Ref. [5] are 7.741 eV/A and 17.874 eV/A, respectively, 
for v = 0.5 a.u., and r s = 1.49 Bohr; and 0.861 eV/A 
and 0.710 eV/A, for v = 0.1 a.u. and r s = 3.04 Bohr. 



